/**
 * @file   test_solver.cpp
 * @author HirasawaYui <yui@Ubuntu18-04>
 * @date   Thu Nov 19 10:32:21 2020
 * 
 * @brief  
 * 
 * 
 */

#include "MultigridSolver.h"
#include "FEMSpace.h"
#include <typeinfo>
#include <math.h>
#include <ctime>
#define pi 4.0*atan(1.0)
double f(double *p)
{
    return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
}

double bc(double *p)
{
    if(p[0] == 0)
    {
        return sin(pi * 0) * sin(pi * p[1]);
    }
    else if(p[0] == 1)
    {
        return sin(pi * 1) * sin(pi * p[1]);
    }
    else if(p[1] == 0)
    {
        return sin(pi * p[0]) * sin(pi * 0);
    }
    else if(p[1] == 1)
    {
        return sin(pi * p[0]) * sin(pi * 1);
    }
}

double f2(double *p)
{
    return 2 * pi*  pi* cos(pi * p[0]) * cos(pi * p[1]);
}

double bc2(double *p)
{
    return cos(pi * p[0]) * cos(pi * p[1]);
}

int main(int argc, char* argv[])
{
    //g++ -o main test_solver.cpp -std=c++11 -I /usr/include/eigen3/ -g
    int n = 6;    
    RectangleDomain* r = new RectangleDomain({{0,0},{0.5,0},{0.5,0.5},{0,0.5}});
    //Mesh<2>* m = new P1Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Triangular_1_Element();
    //Mesh<2>* m = new Q1Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Quadrilateral_1_Element();
    //Mesh<2>* m = new Q2Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Quadrilateral_2_Element();
    Mesh<2>* m = new P2Mesh(r,{POW2(n), POW2(n)});
    Element<2>* e = new Triangular_2_Element();
    Equation<2>* equ = new PossionEquation<2>();
    equ ->SetBoundaryConditionFunction(bc);
    equ ->SetRightHandsTermFunction(f);
    Possion_2D possionproblem(m,e,equ);
    possionproblem.AssembleStiffMatrix();
    possionproblem.AssembleRightHandsTerm();
    //possionproblem.DealWithBoubdaryCondition();
    //P1_MGSolver<Eigen::SparseMatrix<double> > Solver(possionproblem.A(), possionproblem.Rhs(), bc, m, n);
    //Q1_MGSolver<Eigen::SparseMatrix<double> > Solver(possionproblem.A(), possionproblem.Rhs(), bc, m, n);
    //Q2_MGSolver<Eigen::SparseMatrix<double> > Solver(possionproblem.A(), possionproblem.Rhs(), bc, m, n);
    P2_MGSolver<Eigen::SparseMatrix<double> > Solver(possionproblem.A(), possionproblem.Rhs(), bc, m, n);
    Solver.compute(possionproblem.A());
    Solver.setTolerance(1e-12);
    Solver.Solve();
    return 0;
}
    
